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The phase-turbulent (PT) regime for the one dimensional complex Ginzburg-Landau equation 
(CGLE) is carefully studied, in the limit of large systems and long integration times, using an 
efficient new integration scheme. Particular attention is paid to solutions with a non-zero phase 
gradient. For fixed control parameters, solutions with conserved average phase gradient v exist 
only for \v\ less than some upper limit. The transition from phase to defect-turbulence happens 
when this limit becomes zero. A Lyapunov analysis shows that the system becomes less and less 
chaotic for increasing values of the phase gradient. For high values of the phase gradient a family 
, of non-chaotic solutions of the CGLE is found. These solutions consist of spatially periodic or 

• aperiodic waves travelling with constant velocity. They typically have incommensurate velocities 

for phase and amplitude propagation, showing thereby a novel type of quasiperiodic behavior. The 
main features of these travelling wave solutions can be explained through a modified Kuramoto- 
■ Sivashinsky equation that rules the phase dynamics of the CGLE in the PT phase. The latter 

CO ' explains also the behavior of the maximal Lyapunov exponents of chaotic solutions. 

O ' 

o 

oo 
O 



I. INTRODUCTION 

The complex Ginzburg-Landau equation (CGLE) plays a fundamental role in the study of spatially extended 
systems. It describes the dynamics of a generic spatially extended system which undergoes a Hopf bifurcation from 
a stationary to an oscillatory state. Sufficiently close to the bifurcation point, any such system can be described by 
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the CGLE §. For a review of experimental systems described by the CGLE see 

In the present paper we want to discuss only the CGLE in one spatial dimension. Even though this lacks important 
features of the CLGE in higher dimensions, it displays a variety of chaotic and stationary regimes for different values 
of control parameters. Recently, Shraiman et al. || and Chate j| have provided a careful description of these regimes. 
In particular, four different chaotic states have been identified: the phase-turbulent (PT), the defect-turbulent (DT), 
the bichaotic and the spatiotemporal intermittent regime Q. The fundamental states are the PT and the DT regimes. 
The bichaotic state is characterized by the fact that PT or DT chaos can arise, depending on the initial conditions 
— while DT chaos or stable (periodic) solutions can appear in the intermittent region. The PT and DT phases have 
been the object of several detailed studies |§J^^J^|. 
Let us write the one-dimensional CGLE as 

A t = (l + i Cl )A xx + A-(l-ic 3 )\A\ 2 A (1) 

where the parameters c\ and C3 are real positive numbers, while A(x,t) = p(x, t) exp[iip(x, t)] is a complex field of 
amplitude p and phase ip. Eq. (|l]) reduces in the limit (01,03) — > (00,00) to the integrable nonlinear Schrodinger 
equation, an equation well known to admit soliton-like solutions ||. In the opposite limit, (01,03) — > (0,0), the real 
Ginzburg-Landau equation will be recovered, an equation related to symmetry-breaking instabilities of non-oscillatory 
type §• 

Equation (1) admits plane- wave solutions of the form 



Aq(x, t) = VI - k 2 exp[i(kx ~ fl t)} (2) 

where Oq — ~ c 3 + ( c i + c 3)k 2 - Below the so-called Benjamin- Feir (BF) line defined by 03 = 1/ci these solutions are 
linearly stable provided k 2 < (1 — ciC3)/(2(l + c§) + 1 — C1C3) M. Outside this "Benjamin-Feir band", one has an 
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Eckhaus instability against long-wavelength modes. Slightly below the BF line one finds multistability, i.e. the final 
solution depends on the initial conditions ||. Above it, i.e. for C3C1 > 1, all plane wave modes are unstable. The 
latter is the region we are interested in. 

For generic (not necessarily periodic) solutions, the average phase gradient is defined as 

v = — / dxd x ip(x, t) , (3) 
L Jo 

where L is the system size (we assume periodic boundary conditions A(x + L,t) = A(x,t) throughout this paper, 
unless stated differently). The PT phase is defined as that regime where v is conserved [[|. It is encountered just 
above the BF line (i.e. for C3 > 1/ci). It is a state where the chaotic behavior of the field is essentially ruled by the 
phase dynamics. The amplitude, always bounded away from zero for Eq. ([|) to make sense, shows small fluctuations. 

This is different in the DT regime in which the amplitude-dynamics becomes essential ||,^|. Large amplitude 
oscillations are observed in this phase which occasionally drive p(x, t) down to zero. Events where this happens are 
called defects. At a defect, the phase is of course no longer defined, and the appearance of a defect implies that the 
phase difference ip( x + e i t) — ip(x — e, t) jumps by ±27r. Thus v is no longer conserved. 

The transition from the PT to the DT regime happens at a line L\ in the (ci, C3)-plane. To locate this transition 
precisely [|]]7|], one needs some order parameter which allows to distinguish the two phases. Several parameters have 
been proposed recently: the density of defects, the phase and amplitude correlations lengths, and the Kaplan- Yorkc 
dimension. But as observed in R] , the latter are not suitable as order parameters as they do not show any non-analytic 
behavior at the transition line L\. The only good order parameter among the above candidates is the defect density 
Sd- Its value is > in the DT regime, while it vanishes when approaching the PT phase |3|,0]. To present other 
suitable order parameters will be one of the aims of the present paper. 

In the PT regime, all observables depend of course on the value of v, as ergodicity is broken in this phase. Never- 
theless, in literature only few studies have been devoted to solutions with v = [ fH)|JTl] ]. To fill this gap is the second 
goal of the present paper. It will be shown that most observables depend indeed rather strongly and systematically 
on v. 

Since our study is mostly numerical, we need an efficient integration routine. We found that a novel scheme, 
similar to but more efficient than the one introduced in Jl3| ], gave excellent results. This scheme is introduced in the 
next section and compared with the usual pseudo-spectral codes. In Section III the transition from phase to defect 
turbulence is investigated with the help of a new order parameter. A Lyapunov analysis for the solutions with v ^ 
is reported in Section IV together with a detailed characterization of the observed stable solutions. In Section V it 
is shown how a modified Kuramoto-Sivashinsky equation for the phase is able to reproduce the main features of the 
dynamics of the considered solutions of the CGLE. Our results are summerized in Sect. VI. 



II. INTEGRATION SCHEME 



The algorithm applied in this work is a operator splitting scheme similar to the well known "leap frog" algorithm 
for ordinary differential equations (this should not be confused with the "staggered leap frog" |l2| for PDE's). This 
algorithm is based on the fact that the dynamics consists of two independent terms, one nonlinear but local and the 
other nonlocal but linear, each of which by itself can be easily integrated. In the present case, this splitting is chosen 
as 

8 A 

2-=MA + CA (4) 

with 

MA = A-(l-ic 3 )\A\ 2 A , CA= (l + i Cl )|4 • (5) 

ox z 

For a small time step r, the formal solution A(t) = e ( A/ '+ £ ) i A(0) of Eq. (Q) is then approximated by means of the 
Trotter formula 

For a finite time t = nr, n >> 1, we can lump together n — 1 half steps and arrive at 

e (U+C)t ^ e Arr/2 [e Ct e Arr ]n -l e Ct e Arr/2 (?) 
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Notice that we refrain from giving here a formal error estimate. For a chaotic system, the error will of course increase 
exponentially, and the exponent will depend on the maximal Lyapunov exponent A: roughly we expect 6A(t) ~ e cAr f , 
where c is the constant in the 0(t 3 ) correction in Eq. (||). 

Solving the nonlinear evolution e^ T is easy, by considering amplitude and phase separately. We first find 

Pl (x,t + r) = [e^p}(x,t) = [e-^(l/p(x,t) 2 -l) + l]- 1/2 . (8) 

Inserting this into the equation for the phase, the latter can also be solved exactly to obtain 

fafat + r) = [^ r tp]{x,t) = ^(z,0) + c 3 {T + log[p(x, t)/ Pl (x,t + r)]} (9) 

A non-trivial problem appears only for the solution e ct of the linear part since this is non-local. It is only here that our 
algorithm deviates from previous applications of the leap frog algorithm to PDE's Jl4| , p^|JT^jr^ , [T^ | . In these papers, 
the solution of the nonlocal part was performed in Fourier space. This means that two Fourier transformations have 
to be performed for each integration step. It leads to what will be called pseudospectral method in the following. In 
spite of the speed of FFT algorithms, it can be quite time consuming. It can also lead to intolerable round-off errors, 
which was the main reason in |]l3|] to perform this step by integrating with the Greens function of C. More precisely, 
let us denote by K(x,t) the solution of 

dK/dt = CK (10) 
with initial condition K(x,0) = S(x). For the present case it is simply 

K(x, t) = y^e-^ 2 [cos(/3iX 2 ) - i sin^a; 2 )] (11) 
with /3 = /3 r -if3i= [4r(l + ici)] -1 . Then 

[e CT A](x.t) = [ (%K(t,T)A(x-Z,t). (12) 



Since i"T(£, r) can be computed numerically and stored once at the beginning, the algorithm essentially involves one 
convolution for each time step. 

In practice, space has of course to be discretized as well. Let us denote Ai(t) = A(iA, t), where A is the resolution 
of the spatial grid. The most straightforward approach (used indeed in [^3[ ) would seem to start from Eq. (|l2|), and 
replace simply the integral by a finite sum over the central site plus its 2N nearest neighbors, which would lead to 

N 

[e CT A] t {t)^ Y, K(jA,T)A^(t) . (13) 

j=-N 



But this is not the optimal discrete approximation to Eq. (12). To obtain a more precise one, we first notice that since 
e Cr A is a linear operator symmetric under translations and reflections x — > —x, we can write any lattice approximation 
of it as 

N 

[e^AUt) ^Y K ^ A ^ ( 14 ) 

with unknown coefficients Kj. Next we observe that we can decompose A(x,t) into Fourier modes. Let us assume 
that a good approximation is obtained by a superposition of 2N + 1 modes e lkax with k = —N, . . . N, and a some yet 
unspecified positive constant. We then can determine the Kj such that it gives the exact evolution on these modes. 
A straightforward calculation shows that this implies 

K + 2j2 K 3 c os{jAka) = e -(i+^i)tfe 2 a 2 Vfc = 0,...,iV. (15) 
i=i 

Finally, the value of a is chosen such that the L2 norm 

{1-ftTol 2 +2E3XJ \ K i?} 1/2 of K is minimal. With this choice, 
we hopefully minimize the effect of the neglected modes. 
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In the limit N — > oo, the optimal set of Fourier modes are those with a = 2ir/A, and the above scheme should 
give identical results (up to round-off errors) as the pseudospectral method obtained by applying a discrete FFT to 
A t = LA and integrating it in Fourier space. Thus the error committed in our method with finite N cannot be smaller 
than the error committed by the FFT method, and it should tends towards it for N — + oo. 

In order to compare the precision of our algorithm with that of the pseudospectral code, we have evaluated the mean 
square errors accumulated during a fixed total integration time T = 1 for different values of the spatial resolution A 
and of the number N of convolution channels. In all simulations reported below, the time step and the parameter 
of the CGLE have been r = 0.05, c\ = 3.5 and C3 = 0.9. The errors have been estimated by measuring the distance 
(in Li norm) between an orbit obtained by the algorithm to be tested ('test orbit') and a reference orbit obtained by 
the pseudospectral method. In order to reduce statistical fluctuations, the reference orbit was integrated for a time 
>> T, the test orbit was synchronized to it periodically at times t = mT, m = 1,2,..., and the distances building 
up during the intervals mT < t < (m + 1)T were averaged over. 

For estimating the error due to time discretization, we used a reference orbit with the same A as the test orbit, 
but with four times smaller r. To estimate the space discretization error, the reference orbit had the same r, but 
half the value of A (and, in effect, N = 00). Results are shown in Fig. 1. From these data we can draw the following 
conclusions: 

(i) For all studied values of r and A, the time discretization had a much bigger effect than space discretization. We 
could of course have used a smaller value of r (or a larger A) to avoid this, but our choices were motivated by previous 
simulations |§,|]J^]- Nevertheless, we propose that future simulations should use finer time and/or coarser space grids. 

(ii) The convergence of the spatial truncation errors with N is quite fast. Typically, the error is less than twice the 
error of the FFT code for N > 20. This is in contrast to the scheme of jOl, where N had to be chosen larger than 
w 60. Due to the large time discretization error, we could indeed safely use N = 10. 

(iii) The FFT code needs a CPU time which increases logarithmically with the number L/A of mesh points, while our 
algorithm is essentially linear in N. Using the FFT routine C06ECF of the NAG library, the ratio between the CPU 
times needed by our algorithm and the pseudospectral one was 0.24ZV/ ln(L/A). For typical simulations reported in 
this paper (N = 10, L = 1024 — 4096, A = 1/2), this ratio is w 0.25 — 0.3. Thus our algorithm is indeed considerably 
faster than the pseudospectral one. 

III. ORDER PARAMETERS AND THE PT-DT TRANSITION 

To obtain configurations with non-vanishing v, we used several different initial conditions. A first choice was 

A k (t = 0) = e iukA + (r k + ipk) (16) 

with r k and p k random numbers uniformly distributed in [—0.1, 0.1]. Notice that the randomness is needed to break 
spatial translation invariance. Another ansatz was 

A k (t = 0) = e iukA+s "p k {0) (17) 

where p k (0) was uniformly distributed in [0.2, 1.2], and s k S [—0.01, 0.01]. Finally, in a third type of initial conditions 
we introduced also long range correlations in the amplitude, 

p fc (0) = 0.95p k _i(0) +q k , q k e [-0.05,0.05] (18) 

All these initial conditions give consistent final results: in the regime identified as phase turbulent by previous 
authors, v is indeed conserved for sufficiently small initial values. However, for each pair (ci, C3) there exists an upper 
limit vm above which defects arise in the system, leading finally to a phase gradient v < vm- Therefore, defects will 
arise also in the PT regime, provided that the value of v is sufficiently high. The distribution P{p) of the amplitude 
in the statistically stationary state becomes in general wider with increasing v, as shown in Fig. 2. In particular, the 
lower limit p m in of its support decreases with v. But it seems that this decrease in discontinuous, so that p m i n cannot 
be used as a good order parameter. We will come back to this point later. 

The dependence of vm on C3 is plotted in Fig. 3 for c\ = 3.5. This fixed value of c\ was chosen for ease of comparison 
with previous studies. We see an almost linear decrease, and vm reaches zero exactly at the line L\ which marks the 
onset of DT (this was also seen in |l0||). 

Thus vm can be used as order parameters to characterize the transition from phase to amplitude turbulence. From 
the general theory of continuous phase transitions we would in general not expect a linear behavior of the order 
parameter as the transition is approached. Thus we should be suspicious that the linearity seen in Fig. 3 might give 
way to a power law for very small vm- 
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In order to understand this better and to characterize better the transition between phase and defect turbulence 
(always concentrating on the value c\ = 3.5 which we assume to be a generic point on the line L±), we performed 
several simulations with high statistics. For each value of C3 we performed ~ 50 to 70 simulations with v > vm and 
with independent initial conditions, and followed them for > 10, 000 time units (some simulations were followed even 
for > 10 5 time units), unless a defect is found — in which case v is decreased and the simulation is repeated. We 
carefully checked that vm did not depend on the system size, by using systems with L ranging from 1024 to 4096. 
Nevertheless, the present simulations did not show any significant deviation from linearity in the dependence of vm 
on C3, and our best fit is 

v M ^ -(0.74 ± 0.01)c 3 + 0.57 ± 0.01 . (19) 

From this we can estimate the critical value of C3 to be cj = 0.76 ± 0.03. Notice that the linearity of vm agrees with 
a mean field description, since the critical exponent for the order parameter in cases without spontaneous symmetry 
breaking is 1. 

Previous estimates of Cg (for the same value of c\) were all made by approaching the transition from the DT side. The 
defect density Sd was measured in ||[7j by counting how often the total phase difference 5ip(t) = ip(x + L,t) — ip(x, t) 
changed when going from t to t + r. Assuming a scaling law Sd ~ (03 — C3) , the authors of jj) found Cg ~ 0.77 and 
a ~ 2. In contrast, Cg = 0.74 and a ~ 6.8 was found in [Q. Since the latter used much higher statistics (iteration 
times up to 10 7 units, L = 4096 against L — 1024 in ||) one might consider it as more significant. But the very large 
value of a suggests that the data of presumably do not follow a power law, and the extrapolation to Sd = is 
doubtful. Indeed, another extrapolation was also tried in Q which gave Cg = 0.70. Thus our result quoted above is 
in good agreement with previous estimates, in view of the latter's uncertainties. 

To obtain an independent estimate of Cg and of an eventual scaling exponent governing the approach from the DT 
side, we performed simulations (with v = 0) for several values C3 > Cg. In these simulations we did not measure actual 
defects but near defects. More precisely, we made histograms of P(p) similar to the one shown in Fig. 2. A number 
of such histograms are shown in Fig. 4 (actually, the quantity shown in Fig. 4 is P(p)/(2irp), i.e. the 2-dimensional 
density). We see much wider distributions than in Fig. 2, reflecting again the higher degree of chaoticity for v = 0. 
From these histograms we can estimate the probabilities 

w(p) = [ P dxP(x) (20) 
Jo 

to have an amplitude < p. We did this for several values of p, and plotted the results against C3 — Cg in a log-log plot 
(see Fig. 5). We see a decent scaling law if Cg = 0.734 ± 0.007 (which we take as our best estimate of c|), with an 
exponent ~ 6.3. This large exponent might however suggest that the data should not be described by a power law 
at all |0. In that case, the true value of Cg could be considerably smaller and close to the "alternative" value 0.70 
found in Q. Thus the agreement with previous simulations and with Eq. (|l^ ) is as good as can be expected in view 
of these uncertainties. 

In Fig. 4 it is seen that p m in (the lower limit of the support of P{p)) does not approach zero continuously as 
eg — * Cg — e. Instead, as C3 approaches c| from below, it seems that p m in jumps abruptly to 0. On the other hand, 
the 2-dimensional density P(p)/(2irp) seems to be flat at p = 0, for all C3 > Cg. As a consequence, the distribution of 
local phase gradients vi oc should decay as P(v[ oc ) ~ z/ ; ~ (here we assume that V[ oc ~ 1/p for large vi oc ), in qualitative 
agreement with observations in [|| . As pointed out in , the creation of a defect can be similar to a noise-induced 
widening of an attractor. This would naturally lead to a defect density which decreases exponentially with C3. But 
we believe that the phenomenon looks more like an (interior) crisis ]H| . It is known that crises can lead to power 
laws with very large exponents pb]| or even to characteristic time scales which increase faster than any (inverse) 
power of the distance from the critical point |2lj] . But it seems that little is known about crises in spatially extended 
systems. The simplest possibility would be that for each pair (ci,v) there is a critical circle \A\ = p* . Once this 
circle is penetrated, the orbit can come arbitrarily close to the origin, and a defect can develop. Unfortunately, we 
do not see any hint for such a circle in our data. The next simple scenario would be that there is a rotationally and 
translationally invariant set of tori p = p*(x, ip) (each of which individually is not symmetric), such that a defect can 
appear as soon as the solution is inside any one of these tori. We see no way how to test this scenario. 

IV. LYAPUNOV ANALYSIS AND STABLE SOLUTIONS 

To measure the degree of "chaoticity" of the solutions with v ^ in the PT regime, we measured the maximal 
Lyapunov exponent A. We did this for several values of C3 and for several ^-values in the interval [0, I'm]- Our data 
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are shown in Fig. 6, for C3 = 0.35,0.40,0.60, and 0.70. They were obtained from single trajectories spanning between 
1.5 x 10 5 and 6 x 10 5 time units. This time systematic differences were seen between different initial conditions (we 
used again Eqs. ( |l~6| ) and (|l7|)), indicating broken ergodicity even for fixed v. Thus there seem to be several coexisting 
attractors, with slightly different values of A. This complicates the situation of course, and it might be responsible for 
the slight non-monotonicities seen in Fig. 6. Apart from these, the most evident feature seen in Fig. 6 is the strong 
decrease of A with v. Indeed, it seems that in general A — > for v — > vm- This decrease is seen even more clearly if 
averages are taken over several initial conditions, as is done in Fig. 7 for C3 = 0.5. In this figure, A was averaged over 
15 - 35 trajectories for each value of v, with t = 10 5 for each trajectory. 

The progressive ordering of the states with increasing v can be directly seen when plotting snapshots of individual 
solutions. A sequence of such snapshots is shown in figures 8 (amplitudes) and 9 (phases). While the configurations 
seem chaotic for small u, they are more and more regular for larger v, to become finally periodic for v w vj^. While 
it is not surprising that the measured A was zero for Figs. 8e and 9e, it is surprising that the same is also true for 
Figs. 8d and 9d. Roughly, the latter solution can be described as a sequence of periodic "patches" (each similar to the 
solution shown in Figs. 8e and 9e), interrupted by "faults" in which the amplitude is nearly constant and the phase 
decreases linearly. 

We have observed such non-chaotic solutions for all C3 < 0.5 when v — > vm- They can all be expressed in the form 

A(x, t) = h{x - vt)e l(ux ~ ut) (21) 

where h(£) = p (£) e l ^°^ is in general complex and periodic with period L. Notice that also its phase is periodic, so 
that v in Eq. ( |2l|) is the total average phase gradient. We found both solutions which are periodic with period L 
(forced upon them by the periodic boundary condition), and solutions with smaller period (as, e.g. in Fig. 8e and 9e). 
It is natural to assume that the former would lead to spatially aperiodic solutions for L — > 00. The periodic solutions 
can be considered as Bloch waves where the periodicity is however self generated and not imposed by some external 
potential. 

Correspondingly, amplitude and phase can be expressed as 

p(x,t)=p(0 ; i/)(x,t)=ip ^)-wt + vx. (22) 

where £ = x — vt. While the motion of the amplitude is a simple shift with velocity v, the evolution of the phase is 
more complicated. At fixed £, i.e. in a coordinate frame moving with velocity v, the phase increases linearly with 
dtp/dt\^ = vv — m. But there is no fixed phase velocity, i.e. there is no frame in which the phase is strictly constant. 
It fluctuates around a constant value in a frame moving with velocity (v p h) = u)/v, which is thus the average phase 
velocity. Typically, we found that v and (v p h) are not commensurate. Thus solutions of the type of Eq. (|2l]) display 
a novel type of quasiperiodic motion. Similar waves have been found previously in ]22] , |Io| ] . 

While time dependent simulations are needed to verify the stability of these solutions, their existence can be checked 
much more easily. Indeed, they are obtained by solving an ordinary differential equation for h(£), 

(1 + ici)h 66 + [2ii/(l + ici) + v]h s + [1 - v 2 {l + ici) + iu]h - (1 - ic 3 )\h\ 2 h = (23) 

We should notice that for C3 > 0.6 the picture changes slightly: A still decreases with v, but non-chaotic solutions 
are no longer observed. 

Another way to simulate regular solutions with periodic h(£) is to use small systems where L is equal to its period. To 
obtain arbitrary values of v, we have to use in this case twisted boundary conditions p(x+L) = p(x),ijj{x+V) = ijj(x)+9 
with a twist 9 = vL. By means of such simulations we checked again the existence of solutions (j2l|). We checked also 
that there was an upper limit on v which agreed roughly with vm, though its precise value depended on L, as one 
might have expected: depending on its precise value, an integer number of waves will fit into L with more or less ease. 



V. A MODIFIED KURAMOTO-SIVASHINSKY EQUATION 



The behavior of A as a function of v and the origin of the observed non-chaotic solutions can be explained in the 
framework of the modified Kuramoto-Sivashinsky equation (MKSE) derived in 
the main assumptions in deriving the Kuramoto-Sivashinsky equation (KSE) [^[ 

It is well known that just above the BF line the amplitude is nearly constant, and its dynamics is essentially ruled 
by the phase behavior More precisely, the KSE is obtained by making two main assumptions: first, A{x,t) is 
obtained by adding a small perturbation to the spatially constant solution, A(x, t) = e~ lQot + u(x,t); and secondly, 
u(x,t) is a "slaved" function of ip(x,t) and its spatial derivatives. 



|23|]. For completeness, let us recall 
24| and its modification. 
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In deriving the MKSE we just replace the first assumption by making a perturbation around the plane wave solution 
given in Eq. (||). After a tedious and non-trivial calculation p] one finds that the amplitude is given by 



p(x, t) = y/(l- q(x,t))* + (c 3 q(x,t))* (24) 
with q(x,t) = [ci6%ip(x,t) + (d x ip{x,t)) 2 ]/2, and the phase satisfies 

V- + ^ ] ^ xx + n<£ ] (Vx) 2 + ^ X xxx + n£ 2 = o (25) 

where fig = c i c 3 — 1) ^4 = ( c i + c 3)' = c i(l + c i)/2 an d fi£ = 2ci(I + c 2 ). Notice that the latter is just the 
KSE except for the last term. 

In order to verify if Eqs. (^) and (|24|) well approximate the dynamics of the CGLE in the PT regime, two checks 
have been performed. As a first test, we have evaluated p(x, t) by means of Eq. ( |24| ) in the whole PT regime and for 
several values of v. A comparison with the amplitudes obtained directly from the integration always shows satisfactory 
results, see Fig. 10. Notice that this is even true for the extremal values C3 = 0.4 and 0.7, i.e. for two parameter 
values which are near to the BF line and to the line L\. 

As a second check, we have derived from Eq. ( p5| ) expressions for the amplitude velocity v and the frequency us 
associated to the non-chaotic solutions ( pl| ) . We first observe that the ansatz (^lj) leads to 

w = -i>(jc,t)-vtf>t } (£) ; d x( f>(x,t)=^ (0 , (26) 



where primes indicate derivatives with respect to £ = x — vt. We then substitute this into ( pq ) and integrate over all 
x. Due to the periodic boundary condition some terms drop out and others can be simplified, leaving us with 

c = fi + ^ 2) ((^) 2 }-^ 2) («) 2 ). (27) 
where brackets indicate spatial averages, (•) = L^ 1 dx ■ (a similar result was derived in p3| , but with the last 



term neglected). To obtain an expression for v we multiply both sides of ( |25| ) by if>' (C) before averaging, arriving at 

„ _ (rt^M) 3 + 2K^) 2 ] + siPmWo - KVtf) 2 ]) 



(28) 



Numerical tests for several parameter values are shown in Table 1. In all cases the agreement between the measured 
values and the right hand sides of Eqs. (|27]) and (^8|) is very good. 

Having verified that Eq. (^5|) gives a good approximation for the phase dynamics of the CGLE in the whole PT 
regime, we can now relate our numerical findings to several other observations in the literature. This is again based 
on an observation by Sakaguchi |23| who noticed that Eq. (p5|), once ip x is approximated by its average value v in the 
last term, can be rewritten as an equation for s — — ipx " 

S -l - ^2 &xx 2f^2 ^4 &xxxx ~\~ ^Z^xxx • (29) 



with ^3 = v^a\- This is exactly the so-called Kawahara equation 25 2^] which has been studied numerically |25| and 
theoretically |2^ , |27|| . Periodic wave trains formed by pulse- like structures have been found in this equation. It has been 
argued p{| that they are due to the dispersion term cx s m , a conclusion which was confirmed in a careful analysis by 
Chang et al. p7j. These studies confirm also the coexistence of chaotic and non-chaotic solutions (depending on the 
initial conditions). Chaotic attractors dominate for small values of the dispersion constant ^3, while for increasing 
5^3 periodic attractors prevail. Finally, above a threshold value, only non-chaotic solutions are found. Since ^3 oc is, 
our results about the decrease of A with increasing v and the appearance of wave train solutions for v — > v-^i are fully 



consistent with the findings of |27|. Finally, the observation of Chang et al. that stable solutions exist only below a 
certain threshold — which in our case turns out to be v ~ 0.1 for 0.35 < C3 < 0.5 — explains why no regular solutions 
were found above C3 = 0.6: for C3 > 0.6 the value of vm is below 0.1. 



VI. CONCLUDING REMARKS 



In this paper a new integration scheme has been introduced which is faster than the usual pseudo-spectral code 
at comparable accuracy. It has been successfully employed to study the CGLE in the PT regime for large systems 
(L > 1024) and long integration times (t > 10 5 ). Although we have not applied it to other systems than the CGLE, 
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we believe that it should be useful also for the nonlinear Schrodinger equation and for chemical reaction-diffusion 
equations. 

Particular attention was paid to solutions with non-zero average phase gradient v. For fixed control parameters we 
always found a maximal conserved phase gradient vm- It can be used as an order parameter to characterize the PT 
regime and to describe the transition from phase to defect turbulence: at the transition, i>m vanishes with the mean 
field exponent. Another order parameter with non-zero values in the DT regime gave the same transition point but 
another exponent which is more similar to exponents found in previous papers. In contrast, the minimal amplitude 
values reached in the PT regime do not seem to be useful order parameters, as they are discontinuous at the transition. 

For small values of v, the majority of the observed solutions are chaotic. But the maximal Lyapunov exponent 
typically decreases with v, and stable solutions are more frequent than chaotic ones near vm- This can be understood 
by assuming that the phase dynamics is ruled by a modified Kuramoto-Sivashinski equation. 

The regular solutions can be either spatially periodic or aperiodic. In both cases they consist of traveling waves 
whose amplitude and phase velocities are typically incommensurate. These solutions can be considered as generalized 
Bloch waves with self generated (periodic or aperiodic) modulation. For spatially periodic regular solutions these 
results were confirmed by simulations of small systems with twisted boundary conditions (i.e., periodic amplitudes 
but ip(x + L) = ip(x) + 6). 

Duringthe final write-up of this paper, we became aware of a paper by Montagne, Hernandez-Garcia and San 
Miguel 1 10 1 where results similar to those reported in Sections III have been found. The analysis reported in [Tc| ] 



confirms that vm can be used as an order parameter in the whole (c±, C3)-plane of the CGLE. 
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Fig.l: Mean square errors due to space discretization and truncation of the convolution for our algorithm (symbols), 
and due to space discretization for a time splitting pseudo-spectral code (solid lines). The astcriscs show the errors 
due to the time discretization (they are the same for both algorithms). All the data refer to c\ — 3.5, C3 = 0.9, 
t = 0.05 and N — 1024. The results for our algorithm are reported for several values of the number of convolution 
channels N. The simulations were done with double precision, whence the round-off errors are at ~ 10~ 15 . 



Fig. 2: Probability distribution P(p) for C3 = 0.6,ci = 3.5 and for three different values of v. 0.015 (solid line); 0.046 
(dotted line) and 0.098 (dashed line). 

Fig. 3: Maximal phase gradient vm as a function of the parameter C3 (asterisks). The solid line represent a linear fit 
for the vm data. 



Fig. 4: Probability distributions P(p) similar to those shown in Fig. 2, but for v = 0, divided by 27rp, and shown on 
a logarithmic scale. For all curves ci = 3.5. From bottom to top, the values of C3 are 0.74, 0.746, 0.75, 0.752, 0.755, 
0.758, 0.761, 0.765, 0.77, 0.783, 0.8, 0.82, 0.85, and 0.9. Except possibly for the bottom most, all curves are in the 
DT regime. The simulations closest to c% were done with L — 4096 and extended over > 2 x 10 5 time units. 

Fig. 5: Log-log plot showing the probabilities to find an amplitude < p against C3 — 0.734. Each curve corresponds to 
a fixed value of p, with p — m/20 for the m-th curve from the bottom. 

Fig. 6: Maximal Lyapunov exponents A versus the phase gradient v for c 3 = 0.35 (a); C3 = 0.40 (b); C3 = 0.6 (c) and 
c 3 = 0.7 (d). Two different sets of initial conditions are considered: ( |l6| ) (asterisks) and ( |l7| ) (circles). The data have 
been obtained for system sizes L = 1024 and L = 2048 and for integration times t = 150,000 - 600,000. 



Fig. 7: Maximal Lyapunov exponents < A > for C3 = 0.5, plotted against v. Each value has been obtained averaging 
over many initial conditions and over an integration time t = 100, 000 for each trajectory. The considered system size 
is L = 1024. 



Fig. 8: Amplitude p{x) for c 3 = 0.5 and for various values of the phase gradient v. (a); 0.01 (b); 0.04 (c); 0.09 (d) 
and 0.18 (e). 



Fig. 9: Same as Fig. 8, but for the phase </>(x). 

Fig. 10: Amplitudes p from the direct simulation of the CGLE (solid line) and obtained through Eq. ( |2~i| ) (dashed 
line): a) c 3 = 0.7 and v = 0.06; b) c 3 = 0.7 and v = 0: c) c 3 = 0.4 and v = 0.26; d) c 3 = 0.4 and v = 0.03. 

TABLE I. Amplitude velocity v and the frequency u of non-chaotic solutions (see Eq. (|2~t])) for several values of C3 and v. 
The first numbers in columns 3 and 4 were obtained from Eqs. (^) and (J2S|) , while the directly measured values are given in 
parentheses. 
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C3 


V 


to 


V 


0.35 


0.092 


0.317 (0.317) 


0.723 (0.726) 


0.35 


0.184 


0.218 (0.218) 


1.454 (1.458) 


0.35 


0.306 


-0.021 (-0.021) 


2.563 (2.569) 


0.40 


0.092 


0.365 (0.365) 


0.731 (0.735) 


0.40 


0.153 


0.305 (0.305) 


1.236 (1.241) 


0.40 


0.282 


0.074 (0.070) 


2.372 (2.337) 


0.50 


0.123 


0.428 (0.427) 


1.256 (1.241) 


0.50 


0.184 


0.355 (0.341) 


1.592 (1.650) 


0.50 


0.245 


0.229 (0.223) 


1.933 (1.984) 
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